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ABSTRACT 

We present results from and the analysis of data from MAXIPOL, a balloon-borne experiment designed to 
measure the polarization in the Cosmic Microwave Background (CMB). MAXIPOL is the first CMB experi- 
ment to obtain results using a rotating half-wave plate as a rapid polarization modulator. We report results from 
observations of a sky area of 8 deg 2 with 10 -arcmin resolution, providing information up to £ ~ 700. We use a 
maximum-likelihood method to estimate maps of the Q and U Stokes parameters from the demodulated time 
streams, and then both Bayesian and frequentist approaches to compute the EE, EB, and BB power spectra. 
Detailed formalisms of the analyses are given. We give results for the amplitude of the power spectra assum- 
ing different shape functions within the I bins, with and without a prior Cf B = Cf B = 0, and with and without 
inclusion of calibration uncertainty. We show results from systematic tests including differencing of maps, 
analyzing sky areas of different sizes, assessing the influence of leakage from temperature to polarization, and 
quantifying the Gaussianity of the maps. We find no evidence for systematic errors. The Bayesian analysis 
gives weak evidence for an EE signal. The EE power is 55^ fiK 2 at the 68% confidence level for I = 151- 
693. Its likelihood function is asymmetric and skewed positive such that with a uniform prior the probability of 
a positive EE power is 96%. The powers of EB and BB signals at the 68% confidence level are I8+34 /iK 2 and 
—3 l^jQ /iK 2 respectively and thus consistent with zero. The upper limit of the BB-mode at the 95% confidence 
level is 9.5 /iK. Results from the frequentist approach are in agreement within statistical errors. These results 
are consistent with the current concordance ACDM model. 

Subject headings: cosmic microwave background — cosmology: observations — methods: data analysis — 
polarization 



1. INTRODUCTION 

Observations of the Cosmic Microwave Background 
(CMB) have dramatically enhanced our understanding of 
the universe. The recent focus has been on the detection 
of polarization in the CMB because it provides information 
complementary to what can be learned from the tempera- 
ture anisotropy. The discovery and characterization of the 
polarization not only confirms the cosmological interpreta- 
tion of the origin of the temperature anisotropy and large- 
scale structures, but also improves the accuracy with which 
we measure parameters in our cosmological model, such as 
the epoch of reionization. So far det ection of CMB po - 
larization has been reported by DASI dLeitch et al\ 12005), 
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CBI dReadheadefa/.l l2004t ISievers et ali l2005h. CAPMAP 
(Bark ats et al.\\2004). BOOME RANG (iMontroy et al1l2005l) . 
and WMAP dPage et a/Jl2006i) . Here we report results from 
MAXIPOL with particular emphasis on the data analysis 
procedure and cosm ological results. A companion paper 
dJohnson et a/.l l2006) emphasizes the polarimetric instrumen- 
tation and observations. 

MAXIPOL flew from the NASA Columbia Scientific Bal- 
looning Facility in Ft. Sumner, New Mexico in May 2003. A 
region of about 8 square degrees, with Galactic coordinates 
I between 110.69° and 114.98° and b between 38.75° and 
42.49°, was scanned during a 7.6 hour night scan. This re- 
gion was located near the star Beta Ursae Minoris (/3-UMi). 
The beam size was 10 arcminutes. We present data col- 
lected with 12 polarimeters that have a center frequency of 
140 GHz. Polarimetry was implemented by rotating a half- 
wave plate (HWP) at a frequency of 1.86 Hz and analyzing 
the modulating polarization with a stationary grid. W e refer 
the reader to the companion paper dJohnson et al. 2006) and to 
other publications for a more t horough review of MAXIPOL 
and its predecessor MAXIMA dHananv et q/JteOOOtlLee et al\ 
[200ltlStompor et q/.ll200lt Uaffe et a/.ll200UIWu et al.t OQ 1 a; 
Abro e et a/Jl2004h . 

The characteristics of CMB radiation can be described us- 
ing the four Stokes parameters: the intensity /, the linear 
polarization Q and U, and the circular polarization V. The 
anisotropy in / (also called the temperature T) has been well 
measured. The circular component V can only arise from par- 
ity violating physics and is believed to be absent from the 
CMB (though this has not been experimentally verified). The 
Q and U parameters arise during the last-scattering process 
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and have been the focus of recent experimental and theoreti- 
cal work. Their presence is evidence for the standard scenario 
of the last scattering process and their characteristics carry in- 
formation complementary to the temperature anisotropy. Al- 
ternative to the parameters Q and U, one may express the po- 
larization in terms of E and B, which are curl and divergence 
free polarization tensors respectively. For noise-free all-sky 
data Q and U can be converted to and from E and B exactly. 
Otherwise one can only make a statistical conversion between 
the two. Models of cosmological evolution usually predict the 
CMB in the form of power spectra: the auto-correlations TT, 
EE, and BB, and the cross-correlations TE, TB, and EB. In 
this paper we report on measurements of the Q and U Stokes 
parameters, and the corresponding EE, EB and BB polariza- 
tion power spectra. 

We used two alternative statistical approaches to extract 
quantities of interest from the data, the Bayesian and fre- 
quentist approaches. Both have been used successfully in the 
analysis of cosmological data. For a given set of data, the 
Bayesian approach gives the smallest error interval for the es- 
timated quantities, but its application is sometimes computa- 
tionally intractable. 

In our analysis we demodulated the timestreams and sub- 
sequently used a Bayesian maximum-likelihood approach to 
make best-fit maps of the Stokes parameters. We then es- 
timated the EE, BB and EB polarization power spectra of 
the CMB using both Bayesian and frequentist approaches. In 
both cases we accounted for known instrumental effects and 
conducted tests for systematic errors. 

This paper is organized as follows. In Section|2]we describe 
our analysis formalism and procedures for estimating Q and 
U maps from the time-ordered data (TOD), and power spectra 
from these maps. In Section[3] we present our results, includ- 
ing the maps, the power spectra, and systematic tests. Our 
conclusions are given in Section|4] 

2. FORMALISM FOR DATA PROCESSING 
2.1. Time-domain processing 

The time-ordered data (TOD) were flagged for the presence 
of transient signals and calibrated using laboratory data and 
observations of Jupiter. They were cut into segments sep- 
arated by gaps longer than 30 seconds, depending on both 
the flagging of transient signals and the stationarity of noise. 
Various tests were performed to ensure the Gaussianity an d 
stationarity of the noise within each segment (Collins 2006). 
From each segment we estimated and removed an instrumen- 
tal signal that was synchronous with the rotat ion of the HWP , 
which we call HWP synchronous signal (Johnson 2004), 
and de convolved the instrumental filters. See iJohnson et al\ 
( 2006]) for more details about these data processing steps. 

The TOD of each of the polarimeters can be modelled as 
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d, = d(t), t is time, sf = T(x,), s, = Q(x,), sf = U(x,), n, is 
the instrumental noise at time t , x, = x(f) is the sky position 
of the pointing at time t, and e is the modulation efficiency 
of the polarimeter. The units of d, are /iK and thus a calibra- 
tion factor converting from the measured voltage to tempera- 
ture has already been included. We present our results in the 
WMAP convention with the Stokes parameters /, Q, and U, 
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Fig. 1. — HWP rotation angle /3 (upper panel) and scan angle r\ 
(lower panel) from a typical subset of the data. 

taking the North Galactic Pole as the d irection of reference for 
the polarization (Hinshaw et al. 2003). The angle a, = a(x,) 
is the rotation angle of a vector pointing along a great circle 
to the zenith, measured relative to the polarization reference 
vector on the sky. The transmission axis of the polarization 
analyzer is oriented at 90 degrees to the zenith direction. The 
angle (3, = /3(x,) is the rotation angle of the HWP relative to 
the transmission axis of the polarization analyzer. During the 
observations, a, changed at a rate of 15 degrees per hour giv- 
ing a frequency of f a « 1.16 x 10~ 5 Hz while (3, varied at 
fp = 1 .86 Hz. Thus U w Afp = 7.44 Hz. The temporal data d, 
are sampled at intervals of At = 4.8 x 10~ 3 seconds. 

The telescope tracked the guide star /3-UMi while scanning 
in azimuth by 2 degrees peak to peak at a constant frequency 

~ 0.06 Hz for the majority of the data. Here the subscript ?y 
denotes the scan angle. Figure[T]shows typical measurements 
of (3 and r\ from a subset of the data. 

To obtain the Time-Ordered Polarization Data (TOPD) we 
demodulated the TOD to produce independent data streams 
for Q and U. Because of the combination of the HWP rotation 
and the sky scan, the CMB signal is in side-bands of the fourth 
harmonic of j3 t . Multiplication of the appropriate sinusoid and 
applying a band-pass filter gives the TOPD from the TOD 
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The brackets denote a band-pass top-hat filter between 
0.05 Hz and 1.5 Hz. This filter selects the frequencies where 
signals are expected and removes residuals of the HWP syn- 
chronous signal, if they exist. Figure [2] shows a power 
spectrum of a section of the data d, before demodulation. 
The residual peaks (marked by arrows) indicate residuals of 
subtraction of the HWP synchronous signal. The residuals 
are at harmonics of fg . (Figure 6 in the companion pa- 
per |JohnsoneraI] (|2006) shows a similar power spectrum for 
a different section of the data where the subtraction is more 
complete.) Figure [3] shows the power spectrum of the data d, 
multiplied by -2cos(/),/e, but before the band-pass filtering; 
that is, the quantity inside the brackets of equation ([3]). The 
gray area indicates the signal band selected by the band-pass 
filter. The low side of the band-pass is determined by the scan 
speed and the high side by the ~1.3 Hz cut-off due to the 
beam. The power spectrum within the band is consistent with 
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Fig. 2. — Power spectrum for a subset of the TOD d t , after removal 
of the HWP synchronous signal (and after deconvolution of instru- 
mental filters). The solid line is the raw spectrum, and the dashed 
line is its smoothed estimate. The arrows indicate multiples of the 
rotation frequency of the HWP fp showing residuals of the HWP 
synchronous signal. We give a more com plete discussion of this 
synchronous signal in Johnson et al. (2006). 
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Fig. 3. — Power spectrum of the demodulated Q of the data in 
Figure [2] before the band-pass filtering; this is the quantity within 
the brackets of equation (3j. The solid line is the raw spectrum, the 
dashed line is its smoothed estimate, and the shaded area indicates 
the pass band of the filter. The spectrum is consistent with white 
noise within this pass band. 

white noise and residuals of the HWP synchronous signal are 
out of band. Simulations with different plausible power spec- 
tra for the underlying signal show that the band-pass filtering 
reduced the RMS of the Q and U maps by 4% regardless of 
the details of the spectra. This loss of signal was compensated 
for by a Monte-Carlo approach, which we describe later. 

Several effects can bias our estimation of the TOPD when 
using equations d2.lt . The beam convolution inevitably re- 
moves signal on angular scales smaller than ~ 10 arcminutes. 
We rectified this by a deconvolution procedure dur ing CMB 
power spectrum estimation using recipes given in IWu et al] 
(2001b). These recipes also cope with the asymmetry in the 
beams. This deconvolution is discussed in section |2~31 A sec- 
ond effect is that an imperfect implementation of the demod- 
ulation may introduce a bias in the estimation of the CMB 
signal. For example, the band-pass filtering is equivalent to 
a convolution in the time domain so that it induces correla- 
tions in the scan direction while giving rise to some loss of 
signal. A third effect is that the deglitching of the data for 
the transients creates small gaps in the TOD, and thus may 



have influenced the demodulation process. The second and 
third effects were estimated and corr ected by a Monte-Carlo 
approach, as described in section l2~3l 

2.2. Map making 

We employed a standard maximum-likelihood method to 
obtain maps of Q and U from the TOPD. In the time domain 
the TOPD can be modelled as 

df=sf+rf, (5) 
where X = Q or U, sf is the CMB signal and nf is the instru- 
mental noise. We note that nf and nf are independent and 
thus uncorrelated because the demodulation processes to ob- 
tain df and df employ orthogonal kernels (see Eqs. ( 12. j} ). 
We model the CMB signal as 

s?=A tp $, (6) 

and use the Einstein summation convention when appropriate. 
Here A tp is the pointing matrix giving the weight of pixel p in 
observation t , and s x is the CMB signal in the pixel. We took 
the pointing operator A, p to be unity when observing pixel p at 
time t and zero otherwise. That is, we assumed the signal sf to 
be constant within pixel p. This model of pixelization induces 
an extra convolution effect in addition to that from the beam. 
T o deal with these c onvolution effects we followed the recipes 
in lWu et ~ai\ d2001bl) . These provide a way to transfer all these 
convolution effects into a single Be in multipole space, which 
can then be deconvolved when estimating the CMB power 
spectrum Q. 

With this modeling, we can estimate the pixelized maps m x 

from the temporal data df. In the pixel domain we can also 
model m x as a linear sum of the signal and the noise compo- 
nents: 

m K p =^ p +n x , (7) 

where n x is the noise in the pixel domain. Under the assump- 
tions that the noise in the temporal domain is Gaussian and 
that all CMB maps are a-priori equally likely, the maps m x p 
can be estimated by maximizing the likelihood of the signal 
given the data. This gives 

m x =K x p ,A p , t {K x T x d x , (8) 

where A p ', =Aj pl , N x , = (nfnf,) is the time-time noise corre- 
lation matrix and N pp , = (n x n x , T ) is the estimated pixel -pixel 
noise correlation matrix given by 

N x pp ,= [A pt {N x ,r l A,, p ,Y\ (9) 

We apply equations ([8]) and (O to df and df to give the 
maps and m v p , respectively, as well as the noise corre- 
lation matrices N pp , and N ppl . The numerical implemen- 
tation of these equa tions follows the method described in 
Stom por~ef al.\ (1200 ll) . We use square pixels that are 3 arcmin- 
utes on a side. To simplify notation we construct a column 
vector 

- =(»&)■ (10) 

Note that we use the subscript p to denote pixels in the origi- 
nal maps and the subscript q to denote pixels in the simplified 
notation. Similarly we write 

[N Q , \ 
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Fig. 4. — The sky coverage of MAXIPOL, presented as the number 
of 4.8-ms hits per 3-arcminute pixel. The squares are regions of 
1.1, 1.7, and 2.3 degrees across, which we use to estimate the power 
spectra. Pixels with zero hit count but within the square regions have 
zero weighting in the process of power-spectrum estimation. 

Here the off-diagonal blocks are zero in theory because of the 
orthogonality between np and nf . 
After the maps are formed we apply a filtering to them 
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where W qq i is a filtering matrix. The choice of the filter 
depends on the subsequent step in the data analysis. For 
Bayesian power spectrum estimation the maps are left un- 
filtered (or, equivalently, we apply an identity filter). For the 
frequentist approach we apply a noise-weighting filter to cope 
with the anisotropic noise in the pixel domain. Because the hit 
count per pixel decreases towards the edge of the sky patch, 
as can be seen in Figure|4] the S/N ratio is higher at the center. 
To prevent the resulting power spectra from being dominated 
by the low-S/N pixels, we use 
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where U is a column-vector with all entries equal to unity. We 
also tried other forms of the filter used by other authors (e.g. 
iMontrov et al.\ (120051) ). and found that as long as the biasing 
effect of such filtering is accounted for with the Monte-Carlo 
process that we will describe later, the final CMB power spec- 
tra remain unchanged to a 10% level. For map display pur- 
poses we use a Wiener filter 



where S = S, 
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WF . 



= S(S + NT\ (14) 

T ) is the signal-signal correlation ma- 



trix. Because the Wiener filter amplifies expected signals in 
a model-dependent way it induces more bias than other filters 
that rely purely on the noise level measured from the data. We 
thus do not use the Wiener filter for power spectrum estima- 
tion. 

2.3. Power-Spectrum Estimation 

The Q and U can be expanded into spin-2 spherical har- 
monics in the conventional way: 

(Q ± iU)(n) = {4m ± <n) ±2^m(n), (15) 

£m 

where a E hn and af m are the coefficients for the E- and fi-mode 
polarization, respectively, and n is a unit vector directed in the 



direction of observation. The polarization power spectra can 
then be defined as 

' * (16) 
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where Y and Y' are either E or B. To estimate C Y , we use 
two approaches, one Bayesian and one frequentist, which are 
described in Sec . 12 . 3 . II and 12 . 3 51 respectively. In the follow- 
ing we first lay out the formalism which is general to both 
approaches. 

Because the sky coverage of our observation is finite, we do 
not probe independent C YY for each multipole I. Instead, we 
bin the ts and determine a band power within each I bin. In 
addition, to increase the signal to noise ratio we use three bins 
I = 2-150, 151-693, and > 694 such that only the middle bin 
should have signal given the combination of beam size and 
sky area. 

When estimating the band power and presenting the results 
we must specify a model for the shape of the power spectrum 
within each I bin. Given the binned values, we model the 
power spectra as 



C YY ' =V w W t , b c YY \ 



(17) 



and c Y J 



where the subscript b labels an I bin, both C YY 
treated as column vectors, T>w is a square diagonal matrix 
with diagonal elements equal to 

V u =D e , (18) 

and W(b is a matrix defined as 
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1, when^Gb, 
0, when£^&. 



(19) 



Here c Y b Y is defined as the 'band power', and Di is called the 
'shape function'. 

The model for the shape is encoded in Di . We investigate 
the following four cases: 
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Here the Q(acdm) is the power spectrum predicted by the 
concord ance model of the W MAP+ACBAR+BOOMERanG 
result in ISpergel et al\ (120061) . It is a flat ACDM cosmology 
with n b h 2 = 0.022, n m h 2 =0.13, h = 0.74, r = 0.088, n s = 0.95, 
and erg = 0.74. Consequently the power spectra for different 
models of the shape function are related to the estimated band 
power as 

-rr' 



C 



T->(«) 117 YY' 

„YY' 



(24) 



where n = 1,2,3, or 4, and c b Y n) is the band power to be esti- 
mated. 

Prior to the estimation of the CMB power spectra, we also 
need to determine the effects of the beam convolution and 
of the map pixelization, so that these effects can be decon- 
volved durin g the estimation p rocess. We follow the recipes 
described in IWu et ~al\ (l2001bh . The two convolution effects 
are combined into a single transfer function in the multi- 
pole space: 

"l = -S«(beam)^£(pxl)' (25) 
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Fig. 5.— The B e function of each of the 12 MAXIPOL polarime- 
ters (dotted), of the noise-weighted combination ( solid), of the noise - 
weighted combination for MAXIMA2 (dashed) dAbroe et flHl2004l) . 
and of the map pixelization (dot-dashed). 



where Bj (heam) and Bj (pxl) 



are the effective B\ of the beam and 
the pixelization, respectively. The resulting B\ for MAXIPOL 
are presented in Figure[5] The individual B\ (he ., dm) of the MAX- 
IPOL polarimeters are somewhat smaller than those of MAX- 
IMA and thus produce a wider noise-weighted combination. 
The figure also shows that with a pixel size of 3 arc-minutes 
less than 2% of signal power for I < 700 is attenuated by the 
pixelization. 

2.3.1. Bayesian Approach 

A commonly used Bayesian approach for powe r spectrum 
estim ation is the Newton-Raphson algorithm dBond et ~al\ 
119981) . We attempted to use this method to find the MAX- 
IPOL power spectra, but it failed in a variety of ways that we 
believe were due to the low signal-to-noise ratio of our maps. 
The region of parameter space we are exploring is close to 
the boundary where Q = 0, and the allowed solutions include 
negative power values. This region has a non-smooth like- 
lihood that makes the Newton-Raphson method unreliable. 
We note that future B-mode experiments with a low signal 
to noise ratio in their B -mode detec t ion ar e likely to face the 
same problem. See also lAbroe et al.\ d2004l) for high signal-to- 
noise cases where the Newton-Raphson method is also prone 
to failure. Thus we adopted a Markov Chain Monte-Carlo 
(MCMC) approach. 

The MCMC method explores the likelihood space of the 
map, P(Q|map). It generates lists of samples from a param- 
eter space whose distribution is asymptotically the same as 
the posterior distribution of the parameters. This approach 
has a number of advantages: it fully explores the parameter 
space, making no assumptions about the shape of the likeli- 
hood surface, and can be used in any signal-to-noise regime. 
The MCMC method is valid in cases where the shape of the 
posteriors cannot be assumed to take simple forms, such as 
a Gaussian shape. It also has the important disadvantage of 
high computational cost. 

We used Metropolis-Hastings (MH) sampling, which is one 
of the simplest forms of MCMC. It h as been widely used 
in cosmological parameter estimation dSpergel et al.\ 120061: 
MacTavish et al. 200l iLewis and Bridlell2002l). This method 
is also briefly mentioned in iKovac et al.\ (120021) . In our case 
the parameters are binned Q and the random variable is the 
map vector, m, made up of the Q and U pixel values. We take 
the likelihood of the map vector to be Gaussian, with zero 



mean and inter-pixel covariance matrices C = N+S(Ce), where 
N is the noise covariance and S(Cf) the signal correlation 
matri x of the proposed Q (Tegmark a nd de Oliveira-Costa I 
120011) . Thus the likelihood was calculated exactly using 
logL = -i(m r C _1 m+Tr logC). Our MCMC code was de- 
rived from the spectrum solve r MADspec, which is part of 
MADCAP (iBorrill et aZ.ll2006h . The computation of the in- 
verse matrix was the largest computational step and was done 
with a Cholesky decomposition. 

To generate sufficient samples for estimating the parame- 
ters we used two chains, each of approximately 50,000 sam- 
ples. The calculation required 24 hours on 128 processors 
on the Seaborg supercomputer, which belongs to the National 
Energy Research Scientific Computing Center at Lawrence 
Berkeley National Laboratory. We are in the process of opti- 
mizing the method and believe sufficient computational sav- 
ings can be made to make the algorithm scalable to somewhat 
larger datasets. 

Posterior likelihoods for binned Q can easily be calculated 
with the MH algorithm, since the likelihood of a parameter 
value is proportional to its multiplicity in the chain. To per- 
form these calculations w e used the program GetD ist, apart of 
the CosmoMC package dLewis and Bridle 2002), which also 
performs some convergence tests based on derived secondary 
chains. Additional convergence tests based on power spectra 
of par ameter values in the chain, proposed in iDunklev et al.\ 
(120051) . were also performed. In particular, the variance of the 
mean of the two chains was less than 10% of the mean of their 
variances for the parameters of interest. 

The choice of a good proposal density is critical to optimize 
the convergence rat e of an MCM C chain. We followed the 
recommendations in Lewis (2006): We first re-parameterized 
the space by the eigenvectors of the parameter covariance 
matrix to mitigate the effect of highly correlated parameters. 
New proposed jumps were generated along orthonormal ba- 
sis vectors of this new parameter set which were randomly 
rotated every proposals. The length of the jump was a 
Gaussian random variable with mean zero and the variance 
of the appropriately rotated eigenvalue multiplied by a scal- 
ing factor of 2.4. The covariance matrix was estimated with a 
short non-optimized MCMC. 



2.3.2. Frequenti st Approach 

We also used a frequentist approach to estimate Cj F . Fre- 
quentist ap proach was used in the past in analyses of other 
data sets dHivon et q/.ll2002h iMontrov et a/.ll20"05t iPage et al] 
l2006t iBond et q/.ll2003l) . The following sections describe the 
steps in our analysis and Figure [6] shows a flowchart of the 
pipeline. 

Because the sky patch of our observation is only about 4 
degrees across we approximate it as flat to speed up our com- 
putation. We use a Discrete Fourier Transform (DFT) to ap- 
proximate the multipole expansion. With this approximation 
equation ( TTST ) can be reduced and reorganized as 



b - 



-cos 20k — sin2#k 
sin26*k -cos26*k 



(26) 



where a tilde denotes the Fourier coefficient of the corre- 
sponding quantity, k is a Fourier mode, and 9k is the phase 
angle of k. The multipole number can thus be approximated 
as the wave number, I « k = |k|. As a consequence, equa- 
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Fig. 6. — Flowchart of the frequentist approach for estimating the CMB power spectra and their associated errors. The numbers in brackets 
are the equation numbers associated with the operation. The dashed boxes indicate the operations included in the Monte-Carlo simulation. 



tion ( fT6T > reduces to 



tion d27l > as 



C; 



YY' 



■c, 



YY' 



1 

A 



Y Y'* 



(27) 



|k|=* 




(28) 



\k\eb 



where A is the area of the map used in the DFT in steradians, 
and the brackets denote an average over all the wave vectors 
with |k| =k. 

To increase the S/N ratio per I bin in our results, we com- 
bine the Fs into only three bins and estimate their band pow- 
ers c Y b J n) , instead of the Cj y , with a specified shape function 
(see Eqs. ( 12.31 l and (124-b ). This requires a modification of equa- 



where Dt is the shape function 



k ia me anapc luii^Liun u ( convolved by the mul- 
tipole transform of the sky-coverage window. Note that prior 
to the DFT, we apply the filtering of equations ( fl2b and dT3l 
to the maps. 

Due to the finite sky coverage, the flat-sky approximation, 
and the filtering of maps, equation (f28]l is a biased estimator 
of the band powers. This bias is corrected using the Monte- 
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Carlo approach. 

2.3.2.1 The Pseudo-Band Powers 

We apply the DFT of Equations (l26l i and (l28l l to each of 
the square regions selected from the original maps shown 

in Figure |4] to obtain the estimated band powers cj n ^ , 
where a hat denotes an estimator. We refer to these as the 
'pseudo'-band powers, because they are biased and contain 
noise. These pseudo-band powers can be modelled as 



' EE(S+N)\ 
L b(n) 
BB(S+N) 



) 



{bb}{b>b>} 
„EE(N)\ 



rEE 
C b'(n) 



„EB(S+N) . 
-b{n) 



-b(n) \ 
BB(N) , 
\ C b(n) ) 

\+ H bb<-- H bb' 



-b\n) 



+c; 



EB(N) 
b(n) i 



(29) 



(30) 



where 



H 



(») 

{bb} {b'V }■ 



F (n) -pin) 

p(n) p(«) 
- r bt' +r bl' 



(3D 



We call ±Fj£, in equations ( 12.3.2b the 'DFT transfer ma- 
trix' . It accounts for the biasing effect from the DFT approach 
(Eqs. (O, O, (ESJ, and (EHJ). The + F<;f, and are dom- 
inated by the self-coupling and the geometric mixing of E and 
B polarization respectively dChon et a/.ll2004l) . 

With the formalism (12.3.21 i established, our task becomes 
to obtain an unbiased estimator of cH,%. This requires a inver- 
sion of equations (l29l and (l30l . As will be shown, the forms 
of Gui and ±^ ( "' do not need to be estimated individually. 
Instead, we estimate the overall transfer matrices ^j^Hi///} 
and U-H^y --H { bb ^j given a specified shape function D^f and 
the measured Bi , and then compute their inverses. 



Following the frequentist approach, we estimate the noise 
component c Y b ^ N) in equations ( 12.3.2b by using the previ- 
ously estimated pixel-pixel noise correlation matrix N qq i to 
carry out a Monte-Carlo simulation for the noise-only maps. 
A Cholesky decomposition gives 



N, 



■ LL T , 



(33) 



where L is a lower triangular matrix. Then one realization of 
the simulated noise map is obtained by taking 



(34) 



n i=\ Ja =L 8> 



where g is a vector of Gaussian random numbers with mean 
zero and variance one. Finally, applying the DFT approach 
(Eqs. ( fT2l , ( fT3] l, d26l i, and d28l i) to all these noise maps yields 

YY' (N) 

an estimated c fc(fl) . We use 10,000 realizations for the 
Monte-Carlo to obtain our results. 



2.3.2.3 Unbiased Estimator 



We now construct an unbiased estimator for cSL. Taking 



the inverse operation of equation 

zEE 



bin)- 



Gi'iBjW a i 

±H i b t = ±F i b n e > ,GeiB\Wi b ,. (32) 

Here the subscript {bb} denotes the index of b that runs twice, 

and similarly for {££}. The cUl are the underlying CMB sig- and similarly 

nals of the field, and the are the noise components of 

the pseudo-band powers. The Bi, Gi'i, and ±F b f, will be ex- 
plained as follows. 

We call Bi the 'beam transfer function', which has been 
discussed earlier. It accounts for the convolution effects from 
the beam pattern of each of the polarimeters and from the pix- 
elization during the map-making process. 

We call the Gw in equations (12.3.2t the 'time-domain trans- 
fer matrix' . It is induced by the time-domain processing, in- 
cluding the effects from the demodulation and the deglitching 
of the data. Here we use the same Gw for all the EE, EB, 
and BB modes because a signal-only Monte-Carlo simulation 
of 1,000 realizations indicate that it is the same to an accuracy 
of 2%. This is essentially because the form of Gw is domi- 
nated by the band-pass filter in the demodulation process and 
thus behaves simply as a convolution effect of the sky signal. 



'b(n) 



\"{bb}{b'b'} 



gives 
l 



x 



x EE(S+N)_~EEiN)\ 
-b'in) c b'in) \ 

~BB(S+N) _ A BB(N) I ) 



(35) 



-b'in) 



"b'in) I 



-bin) ■ 



( M <n) 
[+ H bb> - 



-H, 



;EBiS+N) _ ~EBiN)\ 
-b'in) L b'in) J 



(36) 



The inversion operation for Hy}b}{b'b'} (+^bb'~ -^w) 
here is feasible only if the underlying matrices are square, i.e., 
only if the numbers of b's and Z?"s are the same. We thus use 
the same binning strategy for b and b' with only three wide 
bands in obtaining our results. 

To estimate H^-,,^^-,, we employ the following end-to- 
end Monte-Carlo simulation. We inject a unit power into 
c bln) f° r one ^ b m ^ at a time. After multiplying the result- 
ing CjZ) obtained from equation ( 124-b with B£ (beam) , we use 
equation ( fTBI l to obtain the signal-only high-resolution maps 
of Q and U, which are then scanned and processed to pro- 
duce mock TOPD. Maps computed from equation (0 using 
the noise matrices measured from the real data are processed 
through equations (fT2l . (fT3l . d26l i. and d28T i to yield the re- 
sulting band powers. These band powers give one column of 
the transfer matrix H^ bb -,r h , b ,y that corresponds to the chosen 
£ bin for input. A Monte-Carlo simulation with 1,000 realiza- 
tions is used to obtain each of the six columns in H)"?, „,,„, 

{bb\{bb'\ 

Finally an inversion of flS,^, and the previously esti- 
mated noise components are used in equations ( 12.3.2b to yield 
unbiased estimates of the band powers cj^ and thus the power 
spectra C^. 



2.3.2.2 Estimation of the Noise Component 



2.3.2.4 Estimation of Error Bars 
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To estimate the error bars of the power spectra, we again 
employ a Monte-Carlo simulation. First, we simulate maps 
that contain both signal and noise using equations d33l and 
( 1341 ). but with the N qq > replaced by 

C m ,=S m ,+N qq ,. (37) 
Here S qq i is the signal-signal correlation matrix based on the 
cU n y The use of the DFT approach and equations (12.3.2| > 
yields the band powers. We compute 10,000 realizations of 
such band powers, obtain the probability distribution for the 
power value within each £ bin b, and calculate the 68% confi- 
dence intervals. 

3. RESULTS 
3.1. Maps and power spectra 

Figure|7]shows the maps for the Q and U Stokes param- 
eters, the maps convolved with a 10-arcmin Gaussian beam, 
and the Wiener-filtered maps (see Eqs. (fT2l and (TBI ). Note 
that the three sets have different color ranges. Some un- 
cleaned systematic errors in the TOD, such as long-baseline 
gain drift, can manifest as non-Gaussian map components, 
such as scan-striping or gradients. It is clear that none of the 
maps show such visible evidence for systematic errors. We 
will conduct more tests for systematic errors in Section l3~3l 

For determining the amplitude of the power spectra we used 
a square region of 2.3 degrees across centered at I = 113°, 
b = 40.25° (see Fig. U and three £ bins: £ = 2-150, 151- 
693, and > 694. Therefore there are in total nine £ bins 
under consideration, three for each of the EE, EB, and BB 
modes. We report results only for the central bin, unless oth- 
erwise stated. This is the only bin likely to contain signal, 
given the beam size and area of the maps. All the Bayesian 
results are reported after marginalization of the joint poste- 
rior likelihood over un-interesting bins. Because there is no 
tractable Bayesian method to account for the loss of power 
due to the band-pass filtering during demodulation of the time 
streams, we rescaled results of Q for the central I bin by a 
factor of 1.06, calculated from an appropriate frequentist ap- 
proach. Bayesian results are quoted as the mode of the like- 
lihood function with 68% intervals of maximum likelihood. 
Frequentist results are quoted as the median of the probabil- 
ity distribution function with 68% intervals about the median. 
We used a shape function D ( p = 1 /[£(£+ 1)] unless otherwise 
noted. 

TableQ]gives the amplitude of the polarization power spec- 
tra using both analysis methods and for different shape func- 

(n) 

tions D) . For ease of direct comparison we present all re- 
sults in £(£+ 1)Q/(2tt) at the bin center £ = 422. The table 
also gives the results after marginalizing over a calibration 
uncertainty that is assumed Gaussian with er ca i = 13%. The 
Bayesian and frequentist approaches give consistent results. 
Results between different shape functions are also consistent 
within statistical uncertainties. 

We note that our Bayesian results are more dependent on 
shape function than the frequentist results, with a variation 
approximately one sigma. Such results h a ve bee n observed 
by other others, such as in M ontroy et al.\ d2005), where the 
variation is as large as four sigma. In our frequentist results 
the effect is less apparent, mainly because of the Monte-Carlo 
bias correction in our pipeline. The origin of this subtle effect 
warrants further investigation. 

Figure|8]gives the Bayesian posterior likelihoods of the EE, 
EB, and BB modes. Two of the likelihoods (EE and BB) ex- 



TABLE 1 
Amplitude of power spectra 



Shape D ( f 


EE 


EB 


BB 


1/1(1+1) 
ACDM 
l/(2£+l) 
1 


55+ 51 

109+ 130 

01+68 

, , 7 +61 
1 1 ;_ 77 


i££ 

z ^-38 
32+52 


-48^ 
-51+ 37 

-20 
^H+34 



Bayesian approach (inc. <T C al) 



l/£(£+l) 


53+ 57 


14+33 
'1--31 


-30+ 34 


ACDM 




27+]f 


_4!+39 


l/(2£+l) 


88+™ 


29+39 


^n+34 
^-31 


1 


108+80 


34+ 53 

J ^6 


-Al^ 46 



Frequentist approach 



1/1(1+1) 


62+ 52 
oz -45 


3+33 
J -32 


26+ 45 
zo -50 


ACDM 


68+« 


s +34 
-31 


21^ 


l/(2£+l) 


7 o+69 
' -43 


23+31 


38^ 


1 


^11 


o+46 
8 -42 


9 1 +66 
Zi -58 



NOTE. — Amplitude of power spectra 1(1 + X)C\ Y> /(2-k) in /xK 2 for a 

wide band I = 151-693 assuming different shape functions D«"'. Errors are 
68% confidence intervals. The middle block with '(inc. <r ca i)' gives results 
including calibration uncertainty <r ca i = 13%. The first row (bold) is the result 
for which we sh ow t he likelihoods in Figures[S]and|9] and is the result that is 
shown in Figure fTol 

hibit a sharp cut-off at the negative end of the parameter axis. 
The cut-off arises when we marginalize the smooth multi- 
dimensional likelihood over the eight other parameters (bins). 
Our likelihood computation method fails (and consequently 
we have zero samples) when any one eigenvalue of the total 
covariance matrix becomes negative. Since multiple bins con- 
tribute to each eigenmode the cut-offs in the different param- 
eters are coherent, and a sharp edge is formed on marginal- 
ization. Because of such correlations between the cut-offs of 
different parameters, integrating over n parameters can lead to 
a cut-off that scales as C"f x . Therefore the sharpness reflects 
the high dimensionality of the marginalized space. It does not 
bias our analysis. 

The posterior likelihood for EE is skewed positive. By in- 
tegrating the likelihood with a uniform prior over both posi- 
tive and negative values we found a 96% probability that EE 
power is positive. This probability value is unchanged after 
inclusion of a 13% calibration uncertainty. 

The posterior likelihoods for EB and BB are consistent with 
no signal. The 95% confidence intervals for the two modes 
are -53 /xK 2 < £(£+ \)Cf B /{2ir) < 81 ^K 2 and -55 /1K 2 < 
£(£+ l)Cf7(27r) < 57 ^K 2 , respectively. To obtain an upper 
limit for the BZ?-mode we removed the negative region of its 
likelihood (see Fig. [SJ and renormalized the rest. We found 
that at the 95% confidence level [£(£+ l)Cf B /(2 7 r)] 1 / 2 < 10.6 
and 9.5 //K with and without the inclusion of calibration un- 
certainty respectively. 

Figure [9] shows Bayesian two-dimensional joint posteri- 
ors. Joint distributions that include the EE mode are skewed 
positive in EE. There are sharp cut-offs of the joint likeli- 
hood surface for the same reason that they occur in the one- 
dimensional likelihoods. This results in straight edges for 
some of the 95% confidence contours. 

We compare the amplitude of the EE power spectrum with 
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Fig. 7. — The maps of MAXIPOL polarization data (top), the maps convolved with a 10-arcmin Gaussian beam (middle), and the Wiener- 
filtered maps (bottom). 



results from other experiments in Figure [10] The right panel 
in the figure is the same posterior likelihood as shown in 
Figure [8] Our result is consistent with the prediction of 
the concordance model, which has a mean value of £(£ + 
1)C ee /(2tt) = 14 pK 2 for our £ bin. This value falls at the 
65% confidence boundary of our likelihood function around 
the mode. The lighter shaded region in the right panel of the 
figure indicates the 68% confidence region of the posterior 
likelihood. The darker shaded region shows the area under 
the likelihood where EE is negative, containing 4% of the to- 
tal area under the curve. 



3.2. Significance of the measured power 

According to standard cosmological models Cf B and C BB 
are predicted to be about one order of magnitude smaller than 
the Cf E and thus undetected by MAXIPOL. We computed 
the Bayesian posterior likelihoods for £(£ + \)Cf E /(2ir) with 
a prior C EB = C BB = 0. The results with and without calibra- 
tion uncertainty are summarized in Table [2] A comparison 
with Table Q] shows that including the priors gives somewhat 
smaller modes and error bars. 

Table [3] summarizes the confidence level at which the hy- 
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TABLE 3 
Probability for cf E > 



1(1+1) C 7 /27i((iK") 

Fig. 8. — One-dimensional marginalized posterior likelihoods of 
1(1 + l)Cj /(2ir). The sharp cut-offs at the negative end reflect the 
high dimensionality of the marginalized space, and do not bias the 
analysis (see text). 



100 



i4 
3. 



3. 




-100 
100 



100 

(u.K 2 ) 



200 



Fig. 



9. — Two-dimensional marginalized posterior likelihoods of 
£(£+l)C r / /(2tt). The inner (thinner) and outer (thicker) contours 
indicate the confidence regions of 68% and 95% respectively. The 
crosses mark the locations of maximum likelihood. 



TABLE 2 



EE POWER SPECTRUM WITH Cf B = Cf B - 



Shape D<f 


Mode 68% 95% 


Bayesian Cf B = C BB = 


l/W+iyj 

1/(2^+1) 


1 9 440 +83 
lz -21 -38 

■M +59 +130 
H1 -38 -71 


Bayesian C EB = C BB = (inc. a ca] ) 


1/W+l)] 
l/(2€+l) 


19 +41 +94 

iz _ 22 _ 3g 

aq +56 +144 
Hy -48 -80 



NOTE. — Amplitude of power spectrum £(£ + l)Cf E /(2tt) in /^K 2 , and 
68% and 95% confidence intervals assuming different shape functions U?' 
with a prior Cf B = C BB = 0. 



Shape Df 


No Prior 


cf> = cf = o 


i/iw+m 


96% 


83% 


ACDM 


94% 




l/(2£+l) 


98% 


92% 


1 


98% 





NOTE. — Proportion of likelihood with positive Cf for different shape 

functions D«"' with no prior, and with a prior C EB = C BB = (not all shape 
functions were considered in this case). 



pothesis Cf E = is rejected for different shape functions, and 
with and without a prior Cf B = Cf B = 0. The procedure for the 
calculation is identical to the one already discussed in Sec- 
tion 13.11 We assume a uniform prior for all values of EE 
and integrate the area below the appropriate likelihood func- 
tion on the positive side. Because the posterior likelihoods 
are all skewed positive most of the confidence levels for a 
positive Cf E are above 90%. These numbers do not depend 
on the magnitude of the calibration uncertainty because the 
calibration uncertainty is a multiplicative factor, which does 
not change the fraction of area under the likelihood for values 
Cf E > 0. 

3.3. Systematic error tests 
3.3.1. Difference maps 

We divided the TOD into two halves in the time domain 
and processed them separately to yield the CMB maps m^ l) 
and m* 2) . Separate noise correlation matrices were computed. 
The difference maps 



(38) 



were constructed and the resulting noise matrices computed. 
We then estimated the polarization power spectra based on 
these Q and U maps, using both Bayesian and frequentist ap- 
proaches. The rows labeled 'time' in Table|4]show the results 
with and without calibration uncertainty. All these results are 
consistent with zero. 

In a similar manner, we combined half of the 12 polarime- 
ters to make one set of maps, and the other half for another 
set, and computed difference maps and the associated noise 
matrices. The row labeled 'polar' in Table |4] shows the re- 
sults. 

Within statistical uncertainties neither the time-domain dif- 
ferencing nor the polarimeter differencing test gives evidence 
for systematic errors. We note that the sizes of the 68% con- 
fidence intervals in Table [4] are on average larger than those 
of the EB and BB modes in Table [T]because the differencing 
process inevitably increases the noise level per pixel. 

3.3.2. Regions of different sizes 

We also investigated the dependence of the frequentist re- 
sults on the size of the square patch chosen for the power spec- 
trum estimation. The square regions of different size that we 
used are indicated by the boxes in Figured The square region 
of width 1 .7° is centered at I = 1 13.23°, b = 40.2°. The square 
region of width 1.1° is centered at / = 113.13°, b = 40.25°. 
The results are summarized in Tableland are consistent with 
the earlier results. There is no significant increase in the er- 
ror bars when using a smaller region of the maps because the 
edges of the square regions are noisier than the central portion 
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Fig. 10. — The £(£ + Y)Cf E /(2tv) at £ < 1000 from all experiments whi ch reported CMB po larizatio n. Our result is tak en fro m the first row 
of Ta ble □ ( numbers in bold typ e). The other results are from C APMAP dBarkats et 0/1120041) . DASI dLeitch et g]1l2005l) . B03 dMontrov eTal\ 
120051) . CBI dSievers et aZll2005l) . and WMAP dPage et g/.l2006h . The right panel is th e posterior likelihood of our result. The solid curve in the 
left panel is a flat ACDM model from the WM AP+ ACB AR+B OOMERanG result in lSpergel et al.\ d2006h (see text for details). It has a mean 
value of £(£ + l)Cf £ /(2ir) = 14 (iK 2 for our bin at £ = 151-693 (shown as an arrow in the right panel). The shaded regions in the right panel 
indicate the 68% confidence region of the posterior likelihood and the region with negative EE power, which is 4% of the total area below the 
likelihood curve. 



TABLE 4 

POWER SPECTRA OF DIFFERENCE MAPS 



Test Appr. EE EB BB 

time B -19^ -55+|| 8+|> 

time B( CTcaI ) -10+J -55^ 8 +49 

time F 21+^ 8+« - 25 +« 

polar F -gff 35$ 30!« 

NOTE. — Amplitude of power spectra £(£ + \)C yy /(2n) estimated from 
difference maps that were constructed with different divisions of the data 
(column 'Test'). Division was done either in time or by polarimeter (la- 
beled 'time' or 'polar' respectively, see text). Results are given both for the 
Bayesian (B) and frequentist (F) approaches. Error values indicate 68% con- 
fidence intervals. 



TABLE 5 

POWER SPECTRA OF DIFFERENT SKY SIZES 



Region Size 

(x° xx°) 


EE 


EB 


BB 


x=\.l 




12 +52 


99+50 


x=\A 


63^1 


8+50 


14+51 
14 -44 



NOTE. — Amplitude of power spectra £(£ + l)C yy /(2n) estimated from 
square regions of different size (see Fig. |4). The error numbers indicate the 
68% confidence intervals. 

(see Fig.|U) and pixels near the edges have negligible statisti- 
cal weight in the power spectrum estimation. 

3.3.3. Gaussianity test for the maps 



Gaussianity in the pixel-domain signal is an essential as- 
sumption for the methods of power spectrum estimation that 
we used. To test our Q and U maps we applied the Kol- 
mogorov test to the eigenval ue-norma l ized K arhunen-Loeve 
coefficients, as performed in IWu et ~al\ d2001al) . If the signal 
is Gaussian, then the K-L coefficients should be normally dis- 
tributed. In the process, we found that some of the eigenval- 
ues of the noise-whitened signal matrix were negative owing 
to the high noise and imperfectly estimated signal in those 
modes. We thus excluded these modes from the test, but in- 
cluded all the other modes. These coefficients passed with 
a clear margin the Kolmogorov test for Gaussianity at 95% 
confidence. 

3.3.4. Beam asymmetry and polarization leakage 

In certain circumstances an asymmetry in the beam may in- 
duce spurious polarization signals. For example, if an asym- 
metric beam rotates simultaneously with the HWP, the result- 
ing EE or BB spectrum will contain power leakage from the 
TT mode. 

Scans of Jupiter were used to quantify leakage from T to 
Q and U. Jupiter has an inherent pol arization of less than 
0.2% at 140 GHz dClemens et al.W 19901) . which is small com- 
pared to the noise on Q and U during beam mapping. Out of 
12 polarimeters only two showed an instrumental polarization 
signal at a level of 4% and 5%. No other polarimeter showed 
leakage from temperature to Q or U at a level larger than about 
1%, which was the typical noise level for this measurement. 

To quantify this effect on the power spectrum we performed 
an end-to-end simulation. Taking the ACDM as the underly- 
ing model for T we conservatively assumed 3% leakage into 
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each of Q and U, which is equivalent to 4.2% instrumental po- 
larization, for all 12 polarimeters. The TOD were constructed 
using the beam patterns as measured in flight. We processed 
this signal-only TOD to obtain maps and power spectra. A 
Monte-Carlo simulation showed that out of 1,000 realizations 
the largest contribution of this leakage into the final EE or BB 
spectrum was 3 /iK 2 for our main bin i = 151-693; the mean 
leakage was smaller. This test demonstrates that the final re- 
sults were not affected by asymmetries in the measured beam 
profiles and by the polarization leakage. 

4. CONCLUSIONS 

We discussed the analysis of CMB data that were taken with 
a HWP polarimeter. Demodulation of the time-domain data 
based on the rotational position of the HWP gave the Q and 
U data. These data showed a white-noise spectrum at fre- 
quencies well below 50 mHz for the majority of the data at 
a level consistent with detector noise. Most of the data were 
also Gaussian and stationary. We made Q and U maps using 
a maximum-likelihood technique. The maps were also shown 
to be Gaussian and there was no visible evidence for system- 
atic errors. 

We calculated EE, EB, and BB power spectra using both 
Bayesian and frequentist techniques. The Bayesian results 
gave weak evidence for EE power that is consistent with 
ACDM cosmology and with previous results. There was no 
detectable signal for the EB and BB spectra. Results from the 
frequentist analysis were consistent with the Bayesian ones. 
We calculated results for different shape functions and with 
different priors and found that the significance of detection 
of EE power was between 83% and 98% with most of the 
results giving a probability larger than 90%. We gave results 
with and without marginalization over calibration uncertainty. 
Inclusion of the calibration uncertainty does not change the 
significance of detection of EE power. 

We presented results from tests for systematic errors in- 



cluding differencing maps, processing sky regions of differ- 
ent sizes, assessing Gaussianity, investigating beam asymme- 
try and searching for polarization leakage. None of the tests 
showed evidence for systematic errors. 

MAXIPOL is the first experiment to produce CMB data us- 
ing a modulating HWP. The techniques we developed to ana- 
lyze such data should have broad applicability for future CMB 
experiments that are planning to use similar modulation tech- 
niques. 
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